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Abstract —Efficiency, comfort, and convenience are three 
major aspects in the design of control systems for residential 
Heating, Ventilation, and Air Conditioning (HVAC) units. In 
this paper we propose an optimization-based algorithm for 
HVAC control that minimizes energy consumption while main¬ 
taining a desired temperature in a room. Our algorithm uses a 
Computer Fluid Dynamics model, mathematically formulated 
using Partial Differential Equations (PDFs), to describe the 
interactions between temperature, pressure, and air flow. Our 
model allows us to naturally formulate problems such as 
controlling the temperature of a small region of Interest within 
a room, or to control the speed of the air flow at the vents, which 
are hard to describe using finite-dimensional Ordinary Partial 
Differential (ODE) models. Our results show that our algorithm 
produces significant energy savings without a decrease in 
comfort. 

1. Introduction 

Heating, Ventilation, and Air Conditioning (HVAC) sys¬ 
tems are complex mechanical devices that control the climate 
of all kinds of buildings, large and small, residential and 
commercial. In most situations, HVAC systems are used 
to maintain comfortable temperatures, while limiting both 
humidity and air speed away from undesirable levels, as 
described in standards such as ASHRAE 55 [1]. Yet, HVAC 
systems are typically controlled in a centralized or static 
fashion, disregarding variations in building configuration 
(windows or doors opened or closed dynamically), human 
activity, or even human perception of environmental condi¬ 
tions. In this paper we present a control framework, based 
on a PDE-constrained optimal control problem, that takes 
into account localized conditions at the room or even person 
scale. We also present simulated scenarios showing that 
our framework produces significant energy savings when 
compared to classical control strategies, and that it can react 
to changing occupation conditions. 

There are many optimization-based studies of the control 
of HVAC units. Goyal and Barooah [2] studied in detail 
the use of RC network circuits to model the temperature 
within buildings. Kelman and Borelli [3], as well as Hazyuk 
et al. [4], [5], used low-dimension ODE-based models to 
control a HVAC unit using Model Predictive Control (MPC). 
Aswani et al. [6] used a learning-based MPC algorithm [7] 
to account for unmodeled dynamics and disturbance in ODE 
models when controlling HVAC units. Domahidi et al. [8] 
and Fux et al. [9] also used a learning-based method and 
MPC, the first using ADABOOST to estimate uncertainties 
and the second using an Extended Kalman Filter. Ma et 
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al. [10] used Stochastic MPC to handle disturbances in the 
control of HVAC units, also using ODE models. These results 
show that optimal control strategies increase the efficiency 
of HVAC systems, yet the use of ODE (i.e., concentrated 
parameter) models means that there is no detailed control 
of the air flow or of the temperature in arbitrary points in a 
room. 

To address these issues, many authors have implemented 
optimization-methods using PDE (i.e., distributed parame¬ 
ter) models. Even though the dynamic behavior of the air 
flow is mathematically complex due to turbulent dynamical 
responses, its response is laminar in larger areas [11], [12], 
thus it can be analyzed using simpler non-turbulent Computer 
Fluid Dynamics (CFD) models. Moreover, the existence and 
smoothness of the solutions of these non-turbulent models 
has been proved under suitable conditions [13], [14], thus 
they can be used together with gradient-based optimization 
algorithms. Ito [15] studied the theoretical optimal control 
of stationary Navier-Stokes equations coupled with the heat 
transfer equation, finding necessary conditions for the exis¬ 
tence of an optimal argument. 

Different variations of CFD models describing heating and 
ventilation situations in buildings can be found in the litera¬ 
ture, such as the papers by Bathe et al. [16], Sinha et al. [11], 
van Schijndel [17], Waring and Siegel [18], and the book by 
Awbi [19]. Yet, these models are not suitable for integration 
with gradient-based optimization algorithms, which require 
the explicit formulation of all the approximating equations 
and their gradients. It is for this reason that we formulated 
our own numerical discretization of the CFD model using 
the Finite Elements Method (FEM) and mixed boundary 
conditions. 

In recent years. Burns et al. [20], [21], [22] studied the 
optimal control of HVAC systems with linearized Navier- 
Stokes models around an arbitrary steady-state solution, and 
then used a Linear Quadratic Regulator (LQR) controller 
to find the optimal solution. Their results show that PDE 
models, usually considered too complex for online numerical 
calculations, can be effectively used for building control. 
Yet the authors’ use of linearized approximations mean 
that the optimal solution is accurate only if it is close to 
the original steady-state linearization point, which poses a 
serious limitation in practical applications. In this paper 
we use a nonlinear Navier-Stokes model, which means that 
our controller can produce large variations in the climate 
variables, at the cost of more expensive numerical compu¬ 
tations. As we show below, a careful choice of parameters 
and numerical algorithms allows for the nonlinear numerical 
computations to be performed within the time constant of 



this dynamical system. 

This paper presents two contributions. First, we use a 
nonlinear non-turbulent Navier-Stokes model together with 
a convection-diffusion heat equation to model the climate in 
a building. We discretize this model using a Finite Element 
Method (FEM), and we use it to find the optimal control 
for the FIVAC system in a building. Second, we show via 
simulations that it is possible to design control objectives 
that are functions of the residents’ locations, which greatly 
increase the efficiency of the FIVAC system. 

The paper is organized as follows: Section |I^ describes 
the CFD model and the formulation of the optimal control 
problem; Section III discusses the discretization of the PDFs 
forming the CED model, and the discretization of the optimal 
control problem; and Section IV gives the results of our 
simulated experiments after we incorporate small changes 
in the actuators of the HVAC unit, such as independent 
control of each vent in a room or control of the angle of 
the air flow in a vent. Our results validate our hypothesis 
that sizable energy savings can be obtained by introducing 
small improvements in the actuation of HVAC units, thanks 
to the accurate CED models that describe the dynamical and 
distributed behavior of the climate variables in a building. 


II. Problem Description 

A commonly missing key feature in many physical climate 
models used to control HVAC systems is the ability to 
capture the real-time spatial variability of the temperature 
and air flow, depending on the floor plan and configuration 
of the building (e.g., open or closed door and windows). For 
this reason, we use a CFD model, which explicitly considers 
temporal and spacial variations, to describe the interactions 
between the temperature, air flow, and pressure. We then 
formulate an optimal control problem where our CFD model 
appears as a constraint, and whose objective function aims to 
minimize the energy consumption of the HVAC system while 
maintaining the temperature constant at a desired reference. 
The inclusion of spatial variables in our description of the 
climate variables will not only improve the accuracy of our 
estimations, but will also allow us to naturally formulate 
richer problems, such as focusing only on an specific region 
in a room, as shown in Section ||V] 

In the remainder of this section we introduce in detail the 
CFD model and the optimal control problem. 


A. CFD Model 

The foundation of our model is the Navier-Stokes equa¬ 
tion, which couples temperature with free flow convection (as 
explained Section 8 in [19], among other references). As 
shown in the literature, atmospheric air can be modeled as 
an incompressible Newtonian fluid when the temperature is 
between —20 [°C] and 100 [°C] [23], [24]. Hence, we can use 
the Navier-Stokes equation for incompressible laminar flows, 
together with the convection-diffusion temperature model for 
fluids. 

Throughout the paper we make two major simplifications 
to the CFD model. First, we assume that the air flow behaves 


as a laminar fluid which reaches a steady-state behavior much 
faster than the temperature in the building. As mentioned in 
Section |Ij both laminar and turbulent flows are present in 
general in a residential building, for example, in the area 
around HVAC vents [19]. However, Sun et al. [25] found 
only minor differences between laminar and turbulent models 
in a geometry similar to ours, while turbulent models are 
significantly more complex than laminar models [26], [27]. 
Hence, we consider a stationary Navier-Stokes equation to 
describe the fluid behavior, and a time-dependent equation to 
describe the temperature behavior. Second, we consider only 
two-dimensional air flows moving parallel to the ground. 
This assumption intuitively makes sense since the air flow in 
the top half of a room can be accurately estimated using a 
two-dimensional model, mostly due to the lack of obstacles 
(such as furniture). These assumptions reduce the accuracy 
of our model to some extent, e.g., van der Poel et al. [28] 
compared 2D and 3D Rayleigh-Benard convection simula¬ 
tions for a cylindrical geometry and showed that differences 
arise for Prandtl constants Pr < 1, while our model’s 
Prandtl constant is Pr = 1.2. Yet, both assumptions allow 
us to significantly simplify the computational complexity of 
our CFD-based control design (measured by the number of 
variables and number of equality constraints of the model), 
which in turn allows us to compute results on the order of 
tens of minutes, as shown in Section |IV] 

Fet H C be the area of interest. We will denote the 
boundary of H by dCt. Fet u: H be the stationary air 

flow velocity , and p: H —>^ K be the stationary air pressure 
in n. Also, given T > 0, let Tg: H x [0, T] —>■ K be 
the temperature in H. Then, using the formulation found 
in [29], the convection-diffusion of temperature in H can 
be described by the following PDF: 

dT 

- Va; • {K{x)V^Te{x,t)) -|-u(a:) • VxTe{x,t) = 

= 9TAx,t), ( 1 ) 

where H x [0, T] —>■ K represents the heat sources in the 
room, K\ fl —>■ K is the thermal diffusivity , Vx ■ = 

is the divergence operator , and Vx = g|-) is the 

gradient operator . 

Similarly, the stationary air flow in H is governed by 
the following set of incompressible Navier-Stokes stationary 
PDFs: 

Axu{x)+(u{x)-'Vx) u{x)+-'VxP{x)+a{x) u{x) = 
Re ' p 

= gu{x)-, and, (2) 

Vj; • u{x) = 0, (3) 

where H represents all the external forces applied 

to the air (such as fans). Re is the Reynolds number (which 
is inversely proportional to the kinematic viscosity ), p is the 
density of the air , a: H —>■ K is the friction constant, u{x) ■ 
Vx = ui{x) —h U 2 {x) is the advection operator , and 

is the Laplacian operator . Since we do not 


















model the vertical dimension of fl, we omit the Boussinesq- 
approximation buoyancy term proportional to Te, which is 
typically included on the right-hand side of Q. We modify 
K and a to model obstacles to heat and air flow in fl, such 
as walls, doors, and windows, as described in [30], [31]. 
In particular, when the point x corresponds to a material 
that blocks air, we choose a{x) ^ u{x), which results in 
u{x) « 0, and when the point x corresponds to air then we 
choose a{x) = 0. 

We divide the boundary of 17 to two outlets of the HVAC 
system, denoted by Fq, one air return inlet, denoted by F^, 
and the exterior walls, denoted by F^,. Thus F^ U Fq U F^, = 
dVt. 

We use a mix of boundary conditions to model the effect 
of the HVAC system in the room, as explained below. Let 
n{x) be the inward-pointing unit vector perpendicular to the 
boundary at x G dfl. Hence, the air flow has the following 
boundary conditions: 

• The HVAC unit’s fan sets the air flow at Fq, hence 
u(x) = Uo n{x) for each x G Fq, where Uo > 0 is the 
HVAC fan speed. 

• The airflow at the inlet is not constrained, hence u(x) 
is free for each x G Fi. 

0 The airflow satisfies a no-transverse condition at the 
exterior walls, hence u(x) ■ n{x) = 0 for each a: € F^,. 
The boundary condition for the temperature is Te{x) = Ta 
for each x G 917, where Ta is the atmospheric temperature . 
We apply a boundary condition for the pressure equation 
only at the inlet, setting p{x) = pa for each x G F^, where 
PA is the atmospheric pressure . 

B. Optimal Control Problem 

To control the temperature in the zone 17^ C 17, we aim 
to minimize the following cost function: 

J {Te{x,t)-T*ydx + Xiv^{t)^dt +X2ul, (4) 

where A 1.2 > 0, T* is the reference temperature set by the 
user, and tf is the time horizon. The heater power v{t) and 
the fan speed Uo are our controlled variables, the former 
appearing as gT^{x,t) = v{t) for each x G Qh C 17, and 
the latter appearing as a boundary condition. We formulate 
a PDE-constrained optimal control problem using the cost 
in together with the CFD model in 0-© and its 
boundary conditions as constraints. We also add inequality 
box constraints for all the controlled variables, so they remain 
within safety limits. 

As explained in Section our experiments introduce 
variations to the cost function in Q depending on the number 
of available actuators. Regardless, the goal of regulating 
the temperature will remain the same throughout all our 
experiments. 

HI. Numerical Implementation 

Our numerical implementation of the PDE-constrained 
optimal control problem described in Section [ri-B| is obtained 
by first using EEM to transform the CED model in 0-0 to 


a set of ODEs as described in Chapters 3 and 4 of [32], and 
then using the consistent approximation technique described 
in Chapter 4 of [33] which transforms optimal control prob¬ 
lems (with ODE constraints) into nonlinear programming 
problems. After those two transformations, we use commer¬ 
cially available numerical solvers to And approximations of 


the desired optimal control, as described in Section IV 


A. FEM Discretization 

Among the many discretization techniques for PDEs, EEM 
stands out for being compatible with complex geometries of 
the domain FI. Intuitively speaking, EEM approximates PDEs 
by dividing the domain into polygons, and then finding a set 
of ODEs for each vertex, and possibly each facet, of each 
polygon. The resulting set of ODEs has the property that 
each ODE is dependent only on its neighbors. 

Before we can formally describe the EEM discretization, 
we need to introduce extra notation. Let iF^(f7,K") be the 
set of functions from FI to M" belonging to L^(f7,K"’), 
whose weak derivative is also in L^(f7,IR"’) [34]. Note 
that iF^(f7,K"), endowed with the dot product (/, p) = 
fo fix) ■ q{x) dec, is a Hilbert space. Similarly, we denote 

Let {Wk}^Z\ be a polygonal partition of FI, i.e., 
Ufci = FI, int(kFfc) n int(fFj) = 0 for each k ^ j, 
and each Wk is a polygon. If {ecfc}^^ is the set of vertices 
in the polygonal partition and is the set of nodal 

points, then we define the test functions C 

iF^(17,M), and iF^(f7,M^), where Ny,Nyj G N 

and Ny < TV^,, with the following properties for each 
fc G {1,..., Ny} and each j G {1,..., 

• ^k, Pk, and ipk are continuous; 

• ^k, Pk, and ifk are nonzero only in the polygons 
containing Xk', and, 

• Ckixk) = i’kixk) = 1, (P 2 j-iiyj) = [J], and 
T2jiyj) = [i]- 

Then, from 0-0, and using Green’s Formulas (see Ap¬ 
pendix C.2 in [35]), we get the following Galerkin identities 
(as described in Chapter 3.6 of [32]): 


/ dT \ 

/ - {i^ix)'F7^Te{-,t),\/^^k) + 

+ {u ■ X7xTei-,t),^k) = {9T„i-,t),^k); (5) 


- ^ {'^xU,XIxPi) + {{u ■ Fix) U,Pj) + {F/xP,Tj) + 

+ = {gu,Pj); and, (6) 

{Fix ■ u,tlJk} = 0, (7) 

for each k G Ny}, j G 2Nyj}, and almost 

every t G [0,7/]. 

Now, given Nt,,,Nu,Np G N, consider the linearly 
independent sets of basis functions C 

iF^(17,K), and C H^{Fl,M.'^). Using these basis 









functions we can project the variables of our CFD model 
into finite-dimensional subspaces, i.e.: 

Te{x,t) = u{x) =J2vu,j0j{x), 

i=i 

Np 

p{x) = ^r]p,j^j{x). ( 8 ) 

f=i 

Applying the representations in ^ to the Galerkin identities 
in (0 results in a set of Ny ODEs with state variables 
■ Similarly, applying the representations to ([^-Q 
results in a set of 2 Nyj nonlinear algebraic equations with 
parameters and Ny linear ones with parameters 

{pp j}^i. All these differential and algebraic equations are, 
in practice, parametrized by constants corresponding to the 
inner products between basis and test functions, as well as 
their gradients. We omit the technical details of the final 
set of equations due to space constraints, and we refer the 
interested reader to Chapter 3 in [32] for more information. 

B. Optimal Control Discretization 

After the FEM discretization, we effectively have an 
Differential Algebraic Equation (DAE) optimal control prob¬ 
lem where 0 contributes Ny ODEs, 0 contributes 2Ny, 
nonlinear equality constraints, and Q contributes Ny linear 
constraints. Several extra equality constraints are added due 
to the boundary conditions of the air flow and the pressure, as 
described in Section III-AI The actual number of constraints 
due to boundary conditions depends on the number of ver¬ 
tices in the polygonal partition {Wk}i.Z‘i over the boundary. 

The consistent approximation of this type of optimal 
control problem is studied in Chapter 4 of [33]. We follow 
the procedure described there, i.e., we hrst normalize the 
problem using the technique described in Chapter 4.1.2 of the 
same book, and then we use the Forward-Euler discretization 
method to transform the ODEs into a sequence of equality 
constraints. Again, we omit the technical details of the final 
equality-constrained nonlinear programming problem due to 
space constraints. 

IV. Experimental Results 

We simulated a two-room apartment with a square area 
of interest (e.g., the area were a resident is located). 
Our goal is to show that using zoned control over the 
area of interest produces a significant improvement over 
controlling the temperature over the whole room. The HVAC 
system consists of two heaters and two forced-air outlets 
with variable-speed fans. Moreover, in order to show the 
efficiency and stability of our zoned control algorithm, we 
simulated 18 different scenarios with different areas of 
interest, distributed uniformly over the apartment. 

A diagram of the apartment is shown in Figures [T] and 
The apartment’s dimensions are 5 x 10 [m^], the width of 
all outlets and inlet is 0.5[m], and each of the areas has 
dimensions 2x2 [m^]. The two heaters are denoted by 
(left) and (right), with dimensions 1x1 [m^]. 



Fig. 2. Air flow of the experiment in Section |IV-A| with target area 12, 
i.e., the whole apartment. Average air speed in the apaifment is 0.11[m/s]. 
Walls are shown in shaded black. Two outlets and one inlet of the HVAC 
system are marked as To, F], and T i respectively. 


The fluid mechanics are governed by the constants Re = 
0.05 and Pr = 1.2, with k(x) = 10“^ and a{x) = 0 
when X G n is located in free air, while k{x) — 10“'* 
and a(x) — 100 when x G D is located on or in a wall. 
The atmospheric pressure is pA = 101.3[kPa], and the 
atmospheric temperature is Ta = 23.83[°C]. We set the 
desired temperature to T* = 24.83[°C], the time horizon 
to tf = 300 [s], and we assume the function is identically 
zero (i.e., no fans are inside the room). As explained in 
Section II-B gT^{x,t) = v{t) for each x G &h, 9 T^{x,t) = 
u'(f) for each x G and g{x,t) = 0 otherwise. The 
parameters in the cost function 0 are Ai = 0.002 and 
A 2 = 0.001. The optimal control problem finds the fan 
speeds Uy and u'y for Eq and E^, respectively, and the heater 
powers v{t) and v'{t) for Qh and respectively. We set 
the fan speed box constraints to [O.l, l] [™], and the heater 
power box constraints to [0,5] [kW]. 

We discretized the area into Npi = 452 elements, and 
the number of total nodes is Ny = 227. We used hrst-order 
Lagrangian elements to define the test and basis functions 
^k, ^k, and '0fc, thus Nt^ = Np = 227. We used 
second-order Lagrange elements to define the test and basis 
functions ipk and ipk, thus 2Nyy = Ny = 1696. More details 
regarding our choice of test and basis functions can be found 
in Chapter 3.3.1 of [36]. The ODE discretization time step 
was chosen as At = 10 [sj. 

We calculated the total energy usage as the sum of the 
heater energy usage, i.e., v(t) + v'(t)dt, and the fan 
energy usage as /g*/p ^jp, ||u(x)|| p(x) dxdf. Our results 
were obtained using a 16-core Xeon E5-2680 computer 
running at 2.7[GHz], with 128[GB] of RAM. We wrote our 
code using Python, the FEM discretization was computed 
using tools from the FEniCS Project [36], and the non¬ 
linear programming problem was numerically solved using 
the SNOPT library [37] interfaced using the Opt Wrapper 
library [38]. The computation time ranged between 15[min] 
and 45 [min] for each experiment. 


A. Nonlinear Navier-Stokes Model 

Using the nonlinear Navier-Stokes model described in 0- 
(j^ we simulated two major scenarios, the first where the 
objective function in 0 uses = O, i.e., the target area 
























4.2 



(a) Temperature distribution at time t f with target 
area Q, i.e., the whole apartment. 


(b) Temperature distribution at time t f with target 
area flz- 


(c) Temperature distribution at time t f with target 
area flz- 


Fig. 1. Results of the experiments in Section |IV-A| Walls are shown in shaded black, and heaters are shown in shaded blue. Values are in [°C] with 
respect to Tj\. 



(a) Average absolute temperature 
error, in [°C], within Qz at time 
tf with respect to T*. 



(b) Ratio of energy usage over aver¬ 
age temperature change within 
Qz, in [Wh/°C]. 


Fig. 3. Results of the experiments in Section [lV-A| and |IV-B| Columns: 
(A) nonlinear model with zoned control, (B) nonlinear model without zoned 
control, (C) linearized model with zoned control. Each column shows the 
median (red line), mean (red box), and first-to-third quartiles (blue box). 


is the whole apartment, and the second where we use a 
zoned approach with a smaller which moves around 
the apartment to 18 different locations, as explained above. 
Figure la shows the temperature distribution for the first 
scenario, while Figures [T^ and [T^ show the temperature 
distribution for two of the 18 zoned simulations. Also, 
Figure shows the stationary airflow for the first scenario. 

In the first scenario, where = SI, the optimal average 
absolute temperature error in the apartment was 0.502[°C] 
at time t/, and the ratio of energy usage over average 
temperature change within SI was 1009.2[Wh/°C]. We cal¬ 
culated the same statistics for the 18 different zones S7^, 


which are summarized in Column (A) of Figures 3 a and 3b 


To make both scenarios comparable, we recalculated these 
statistics for the first scenario, this time considering the 
average temperature changes in Sl^ instead of SI, which are 
summarized in Column (B) of Figures and Those 
figures clearly show that using zoned control is significantly 
more accurate and more efficient than heating the whole 
apartment. It is worth noting that the zoned approach requires 
roughly half the energy to change the average temperature 
by 1[°C] in when compared to the first scenario. 

The results in Figure indicate that when the resident 
is near one heater, say 0/^, our algorithm automatically 
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Fig. 4. Absolute en'or in temperature distribution, with respect to FEniCS 
simulation, for the linearized model in Section |IV-B| at time t f with target 
area Qz. Values are in [°C]. 


shut down the other heater, say Qh', as intuitively expected. 
Therefore, if it is possible to localize a resident within an 
apartment, e.g., via Bluetooth beacons or using a sensor net¬ 
work, then we can increase the efficiency of the HVAC unit 
significantly without major modifications to the mechanical 
ventilation system. 


B. Linearized Navier-Stokes Model 

We also computed the optimal control using the linearized 
Navier-Stokes model described in [21], [22]. The linearized 
model has clear advantages over our nonlinear model, in¬ 
cluding a larger set of theoretical results supporting it, and 
a faster computation time. On the other hand, linearized 
models perform well only when the values produced by the 
model are close to the stationary linearization point. 

We ran the same experiments as in the second scenario of 
Section IV-A i.e., controlling the temperature in 18 different 
zones. The statistics for average absolute temperature error 
in Liz, and ratio of energy usage over average temperature 
in Liz, are shown in Column (C) of Figures and 3b Even 
though the energy efficiency is comparable when we use 
linearized or nonlinear models, the accuracy is significantly 
different, with the nonlinear model consistently performing 
better than the linearized model. We believe the difference 
is due to the lack of accuracy of the linearized model. As 
exemplified in Figure there is a large error in temperature 
distribution between the linearized model and an accurate 
benchmark simulation using the FEniCS solver, which is 
consistent with the highly nonlinear behavior of the Navier- 
Stokes equation. 









































































V. Conclusion 

Our results open the door to a large number of exiting 
opportunities to improve the energy efficiency of buildings. 
By making small improvements to existing HVAC units it 
is possible to dramatically increase the efficiency of HVAC 
units without a decrease in human comfort. It is worth noting 
that our results do not require, in principle, the use of ex¬ 
pensive variable-speed fans or variable-power heaters, since 
those control signals can be implemented using switched 
strategies [39], [40]. More importantly, our simplifications 
have allowed us to obtain results in tens of minutes while still 
capturing the distributed behavior of the climate variables, 
which is much closer to real-time applications than previous 
results [24], achieving a good trade-off when compared to 
less accurate linearized Navier-Stokes models. 
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